home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlamc5.f < prev    next >
Text File  |  1996-07-19  |  5KB  |  163 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
  5. *
  6. *  -- LAPACK auxiliary routine (version 2.0) --
  7. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  8. *     Courant Institute, Argonne National Lab, and Rice University
  9. *     October 31, 1992
  10. *
  11. *     .. Scalar Arguments ..
  12.       LOGICAL            IEEE
  13.       INTEGER            BETA, EMAX, EMIN, P
  14.       DOUBLE PRECISION   RMAX
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DLAMC5 attempts to compute RMAX, the largest machine floating-point
  21. *  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
  22. *  approximately to a power of 2.  It will fail on machines where this
  23. *  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
  24. *  EMAX = 28718).  It will also fail if the value supplied for EMIN is
  25. *  too large (i.e. too close to zero), probably with overflow.
  26. *
  27. *  Arguments
  28. *  =========
  29. *
  30. *  BETA    (input) INTEGER
  31. *          The base of floating-point arithmetic.
  32. *
  33. *  P       (input) INTEGER
  34. *          The number of base BETA digits in the mantissa of a
  35. *          floating-point value.
  36. *
  37. *  EMIN    (input) INTEGER
  38. *          The minimum exponent before (gradual) underflow.
  39. *
  40. *  IEEE    (input) LOGICAL
  41. *          A logical flag specifying whether or not the arithmetic
  42. *          system is thought to comply with the IEEE standard.
  43. *
  44. *  EMAX    (output) INTEGER
  45. *          The largest exponent before overflow
  46. *
  47. *  RMAX    (output) DOUBLE PRECISION
  48. *          The largest machine floating-point number.
  49. *
  50. * =====================================================================
  51. *
  52. *     .. Parameters ..
  53.       DOUBLE PRECISION   ZERO, ONE
  54.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  55. *     ..
  56. *     .. Local Scalars ..
  57.       INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
  58.       DOUBLE PRECISION   OLDY, RECBAS, Y, Z
  59. *     ..
  60. *     .. External Functions ..
  61.       DOUBLE PRECISION   DLAMC3
  62.       EXTERNAL           DLAMC3
  63. *     ..
  64. *     .. Intrinsic Functions ..
  65.       INTRINSIC          MOD
  66. *     ..
  67. *     .. Executable Statements ..
  68. *
  69. *     First compute LEXP and UEXP, two powers of 2 that bound
  70. *     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
  71. *     approximately to the bound that is closest to abs(EMIN).
  72. *     (EMAX is the exponent of the required number RMAX).
  73. *
  74.       LEXP = 1
  75.       EXBITS = 1
  76.    10 CONTINUE
  77.       TRY = LEXP*2
  78.       IF( TRY.LE.( -EMIN ) ) THEN
  79.          LEXP = TRY
  80.          EXBITS = EXBITS + 1
  81.          GO TO 10
  82.       END IF
  83.       IF( LEXP.EQ.-EMIN ) THEN
  84.          UEXP = LEXP
  85.       ELSE
  86.          UEXP = TRY
  87.          EXBITS = EXBITS + 1
  88.       END IF
  89. *
  90. *     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
  91. *     than or equal to EMIN. EXBITS is the number of bits needed to
  92. *     store the exponent.
  93. *
  94.       IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
  95.          EXPSUM = 2*LEXP
  96.       ELSE
  97.          EXPSUM = 2*UEXP
  98.       END IF
  99. *
  100. *     EXPSUM is the exponent range, approximately equal to
  101. *     EMAX - EMIN + 1 .
  102. *
  103.       EMAX = EXPSUM + EMIN - 1
  104.       NBITS = 1 + EXBITS + P
  105. *
  106. *     NBITS is the total number of bits needed to store a
  107. *     floating-point number.
  108. *
  109.       IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
  110. *
  111. *        Either there are an odd number of bits used to store a
  112. *        floating-point number, which is unlikely, or some bits are
  113. *        not used in the representation of numbers, which is possible,
  114. *        (e.g. Cray machines) or the mantissa has an implicit bit,
  115. *        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
  116. *        most likely. We have to assume the last alternative.
  117. *        If this is true, then we need to reduce EMAX by one because
  118. *        there must be some way of representing zero in an implicit-bit
  119. *        system. On machines like Cray, we are reducing EMAX by one
  120. *        unnecessarily.
  121. *
  122.          EMAX = EMAX - 1
  123.       END IF
  124. *
  125.       IF( IEEE ) THEN
  126. *
  127. *        Assume we are on an IEEE machine which reserves one exponent
  128. *        for infinity and NaN.
  129. *
  130.          EMAX = EMAX - 1
  131.       END IF
  132. *
  133. *     Now create RMAX, the largest machine number, which should
  134. *     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
  135. *
  136. *     First compute 1.0 - BETA**(-P), being careful that the
  137. *     result is less than 1.0 .
  138. *
  139.       RECBAS = ONE / BETA
  140.       Z = BETA - ONE
  141.       Y = ZERO
  142.       DO 20 I = 1, P
  143.          Z = Z*RECBAS
  144.          IF( Y.LT.ONE )
  145.      $      OLDY = Y
  146.          Y = DLAMC3( Y, Z )
  147.    20 CONTINUE
  148.       IF( Y.GE.ONE )
  149.      $   Y = OLDY
  150. *
  151. *     Now multiply by BETA**EMAX to get RMAX.
  152. *
  153.       DO 30 I = 1, EMAX
  154.          Y = DLAMC3( Y*BETA, ZERO )
  155.    30 CONTINUE
  156. *
  157.       RMAX = Y
  158.       RETURN
  159. *
  160. *     End of DLAMC5
  161. *
  162.       END
  163.